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We carry out a comparative study of electronic properties of 2D electron gas (2DEG) in a magnetic 
field of an infinitesimally thin solenoid with relativistic dispersion as in graphene and quadratic 
dispersion as in semiconducting heterostructures. The problem of ambiguity of the zero mode 
solutions of the Dirac equation is treated by considering of a finite radius flux tube which allows 
to select unique solutions associated with each K point of graphene's Brillouin zone. Then this 
radius is allowed to go to zero. On the base of the obtained in this case analytical solutions in the 
Aharonov-Bohm potential the local and total density of states (DOS) are calculated. It is shown 
that in the case of graphene there is an excess of LDOS near the vortex, while in 2DEG the LDOS 
is depleted. This results in excess of the induced by the vortex DOS in graphene and in its depletion 
in 2DEG. We discuss the application of the results for the local density of states for the scanning 
tunneling spectroscopy done on graphene. 

PACS numbers: 03.65.-w, 73.20.At, 72.10.Fk 

I. INTRODUCTION 

Physics of graphene perception begins when one compares Landau levels in two dimensional Schrodinger and Dirac 
theories. Such spectacular phenomenon as unconventional quantum Hall effect 1 '- is caused by the anomaly of the 
lowest Landau level (LLL)2 which for Dirac fermions in graphene is field independent and can accommodate only half 
the usual number of the states from the conduction band and takes the other half from the valence band. The easiest 
way to accomplish this peculiar feature of the LLL is to solve a pair of the Dirac equations that describe excitations 
near two inequivalent K points of graphene's Brillouin zone. Normally this is done in a constant homogeneous 
magnetic field, although this property of the LLL for Dirac fermions is topologically protected for inhomogeneous 
field configurations and in the presence of ripples^. 

The simplest inhomogeneous field configuration which contains nontrivial Aharonov-Bohm physics can be created 
by an infinitesimally thin solenoid. In practice such magnetic field configuration may be obtained when a type-II 
superconductor is placed on top of graphene or semicondunducting hetero junction hosting a 2D electron gas (2DEG) 
with quadratic dispersion. While graphene devices still have to be fabricated, devices like this with a superconducting 
film grown on top of a semiconducting heterojunction (such as GaAs/AlGa As) hosting a 2DEG have in fact been 
fabricated twenty years ago^ and theoretically well studied (see e.g. Refs. p7Hllh . 

Theoretically a problem of the Dirac fermions in the field of Aharonov-Bohm flux was encountered in the context 
of cosmic strings by Gerbert and Jackiw-i^. While for the solutions of the Dirac equation with nonzero angular 
momentum^ the square integrability requirement specifies which of the two independent solutions should be taken, 
they noticed this is not the case for the zero angular momentum. For the zero momentum there is an ambiguity 
as both solutions are square integrable, but divergent as l/y/r at the origin, where r is the space coordinate. The 
ambiguity of the solution selection is caused by the singular nature of the infinitesimally thin solenoid vector potential 
at the origin. This problem has initiated a vast theoretical literature which addresses interesting aspects related to 
the rigorous treatment of the solutions of the Dirac equation with the Aharonov-Bohm potential (see e.g. Refs. [pH - HH 
for a review). In the condensed matter context the Dirac fermions in the field of solenoid emerged during the study 
of the Dirac-Bogolyubov-de Gennes quasiparticles in the vortex state of d-wave superconductors^ (see alsoi 7 . for a 
review). Due to the divergence of the zero mode solutions theory predicts a formation of nonzero local density of states 
(LDOS) near the vortex center. However, this theoretical prediction based on the Dirac nature of quasiparticles in 
(i-wave superconductors does not agree with the results of scanning tunneling spectroscopy (STS) measurements^ in 
high-temperature superconductors. Finally we mention a related problem of the description of topological defects in 
graphene based on the Dirac equation with a pointlike pscudomagnetic vortex which has also been studied intensively, 
see e.g. Refs. HIM 

The purpose of the present paper is to study the electronic excitations in graphene in the field of the Aharonov-Bohm 
flux and compare them with the corresponding results for 2DEG with a quadratic dispersion. We rely on the existing 
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studies of the Dirac fermions in the Aharonov-Bohm potential, but focus on the specific features of graphene such as 
the presence of two inequivalcnt K points which implies that one should consider the solutions for both incquivalcnt 
irreducible representations of the Dirac 2x2 matrices. Also to avoid unnecessary formal complications we consider 
the physical regularization of the problem modeling a finite radius flux tube created by the Abrikosov vortex. We 
utilize the simplest case of magnetic field concentrated in a thin cylindrical shell of small but finite radius R when 
R — > 21 i 22 . In contrast to high-temperature superconductors Dirac description of the quasiparticles in graphene is 
proven valid under the different conditions. In particular, STS measurements of graphene flakes on graphite 2 ^ exhibit 
the structural and electronic properties expected of pristine graphene such as the development of a single sequence 
of pronounced Landau level peaks corresponding to massless Dirac fermions in a homogeneous magnetic field. We 
propose to perform STS measurements for graphene penetrated by vortices from a type-II superconductor, because 
the Dirac theory predicts rather peculiar behavior of LDOS not expected for the 2DEG with a quadratic dispersion 
of carriers. 

The paper is organized as follows. In Sec. |TT]we introduce the model Hamiltonians and discuss the regularization 
of the Aharonov-Bohm potential used in this work. Sec. Illll is devoted to the nonrelativistic case, and the relativistic 
case is discussed in detail in Sec. IIVI The structure of both sections is the same: we consider the solution of the 
corresponding Schrodinger or Dirac equation, construct the Green's function (GF) with coinciding arguments, obtain 
the DOS and study the behavior of the LDOS. In Sec. [V]our final results are summarized. 



II. MODELS AND MAIN NOTATIONS 



We consider nonrelativistic and relativistic Hamiltonians. The 2D nonrelativistic (Schrodinger) Hamiltonian has 
the standard form 

H s = -^{Di + Dl), (2.1) 

where Dj = Vj + ie/hcAj, j — 1,2 with the vector potential A, Planck's constant h and the velocity of light c 
describes a spinless particle with a mass M and charge — e < 0. The Dirac quasiparticle in graphene is described by 
the Hamiltonian 

H D = -ihv F f3( 11 D 1 + l2 D 2 ) + A/3, (2.2) 
where the matrices /3 and Pjj are defined in terms of the Pauli matrices as 

/3 = (7 3 , 07i = (<n,Cffa). (2.3) 

Here C = ±1 labels two unitary inequivalent representations of 2 x 2 gamma matrices, so that one considers a pair of 
Dirac equations corresponding to two inequivalent K± points of graphene's Brillouin zone. In Eq. (|2.2j) vp w 10 6 m/s 
is the Fermi velocity and A is the Dirac mass (gap) , which is introduced in the Hamiltonian for generality. Note that 
we consider the simplest case when the gap has the same sign for ( — ± 1 [see Ref . for a discussion of more general 
cases]. While tight binding calculations show that the quasiparticle excitations in graphene have a linear dispersion 
at low energies^ and are described by the massless Dirac equation with A = 02&, recent STS measurements revealed 
a mass gap near the Dirac point in a single layer graphene sample suspended above a graphite substrate 2 ^. Since this 
gap and its origin are intensively studied both theoretically and experimentally in the last few years, here we consider 
a generic case with a finite value of A. 

The vector potential of a vortex at the origin directed in the e z direction is 

A(r)= 2^ (rXG2) ' (2 - 4) 

where $ = t]^>q is the flux of the vortex expressed via magnetic flux quantum of the electron $o — hc/e with 
r] 6 [0, l[r^ where the value -q = 1/2 corresponds to the flux created by the Abrikosov vortex. The magnetic field is 
then 

B(r)=Vx A = e^$ £ 2 (r). (2.5) 
The essential difference between the Schrodinger 

H s ip = Eij>, (2.6) 
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and Dirac 

H D ^ = E^ (2.7) 
equations in this case can be seen if one squares the latter: 

- H 2 v 2 F (Df + D 2 2 + iCa 3 [D u D 2 ]) = (E 2 - A 2 )*, (2.8) 



where the commutator 



i[D 1 ,D 2 } = -^-B z (r) (2.9) 
he 



which introduces pseudo-Zeeman term which is related to the sublattice rather than the spin degree of freedom. It 
should be mentioned that in the case of graphene the components of the spinor ^(r) are associated with a sublattice 
rather than a spin degree of freedom. Since the Hamiltonian (|2.2[) originates from a nonrelativistic many-body 
theory, the Zeeman interaction term has to be explicitly added to this Hamiltonian. This resembles the case of the 
nonrelativistic Hamiltonian (|2.1[) which becomes Pauli one when the interaction between the magnetic moment of 
the spin and an external magnetic field is added. In the present paper we do not include the spin degree of freedom 
neither in (|2~T1) nor in (j2~2"j) . 

Eqs. (|2.8p and (|2.9p identify the origin of complication s 22 ^ 7 in the problem with a singular vortex (|2.4[) when a 
singularity in B z {v) occurs at a singular point of the differential equation (|2.8|) . To avoid these complications one 
can consider a vortex with a finite radius flux tub o 9 ' 21 i 22 , i.e. with the magnetic field and vector potential written in 
cylindric coordinates r = (r, <p, z): 

B (r) = ^h(r)e z , A(r) = f^e„, (2.10) 
2n 2n r 

where h(r) is a profile function with a compact support satisfying the normalization J °° drrh(r) = 1 and connected 
to the profile function a(r) by the relation h(r) = a'(r)/r. 

The simplest choice of the field distribution h(r) which regularizes the problem with the solutions solely expressed in 
terms of Bessel functions is a magnetic field concentrated on the surface of the cylinder of radius R, h(r) = 8(r — R)/R. 
Then, the corresponding profile function 

a{r) = 6{r-R). (2.11) 

In the limit R — > we recover the Aharonov-Bohm potential (12. 4[) but avoiding formal complications. As shown in 
Ref. 22 , there is no dependence on the detailed form of h(r) in the limit R —> provided that lim,.^o | () ' h(r')r'dr' = 0. 
In the present paper we restrict ourselves by considering the profile function (|2.1ip . 

III. NONRELATIVISTIC CASE 

In this section we recapitulate the results of Refs. [9lfl0l for nonrelativistic case. They are important not only for 
comparison with relativistic case, but also because the relativistic result is constructed using the nonrelativistic one. 

A. Solutions of the Schrodinger equation in Aharonov-Bohm potential and general representation for LDOS 

In the limit R — ¥ the admissible solution of the Schrodinger equation is always a regular solution which in polar 
coordinates r = (r, ip) takes the form 



^ m (r,<p) = \l—e m *J lm+1ll (kr), meZ, (3.1) 



where J| m +,,|(fcr) is the Bessel function with the wave vector k which is related to the quasiparticle energy E{k) via 
E{k) = h 2 k 2 /2M. 

The eigenfunction expansion for the retarded Schrodinger GF reads 

u im — — rv\ x * 
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or after substituting the wave function (|3.1[) it becomes 

M f°° kdk °° 
C%{T,T>,E + iD) = ^ 2 + zQ E e^-^J lm+nl (kr)J [m+v[ (kr'), (3.3) 

m=— oo 

where g 2 = 2ME/H . Since an analytic continuation of the GF (|3 . 3[) on the imaginary axis in the complex momentum 
plane, q —¥ z = iQ is free of singularities, it is convenient to work with the corresponding GF 

M 

G s n (r,r,Q) = —g v (r,Q), (3.4) 

where 

r°° kdk ^ 

9v{*,Q) = -l E J lW^)- ( 3 - 5 ) 

m— — Co 

In Eq. (|3.4p we already set two arguments coinciding, r = r', because in the present work we consider the DOS only. 
As we will see below, the function g v (r, Q) is also used in the representation of the DOS for the Dirac fermions. After 
the calculation of the GF (|3.5|) is done, the LDOS per spin projection can be found by returning back to the real 
momentum axis 

N^(r,E) = -hmG s ri (T,T,Q^-iq + 0), E =^- (3-6) 

The GF G^(r, r', Q) was calculated in Ref. [H using the contour integration technique. A weak point of this cal- 
culation was discussed in 9 -, where the same method was applied to obtain the ry-dependent contribution to the GF, 
AG^(r, r, Q) = Gf.(r, r, Q) — Gq (r, r, Q) with coinciding arguments when the approach is valid. Referring to the 
derivation of2£, here we simply start from the corresponding expression for Ag v (r, Q) — g v (r, Q) — go{r, Q) obtained 
inRef.[H 

p oo poo r)(v — uj) 

Notice that Eq. (13.7[) coincides with the corresponding expression from 9 up to a coefficient. Changing the variables 
to x — (v — w)/2, y — (v + ui)/2 one can obtain from (|3.7j) the final expression 

r dy r ^H{2y-i)x] __ 2ar 

lo Jo 



A 5 „(r,G) = ^^ / dy dx — e -2Qrcosh3;coshiy (3g) 

' ' coshx 



which we will use in what follows. It turns out that the investigation of the full DOS is simpler than the analysis of 
the LDOS. Thus in the next Sec. lHIBI we firstly consider the DOS and return to the LDOS ([3~6]) below in Sec. IIII CI 

B. The density of states 

The full DOS per spin projection is obtained from the LDOS (|3.6I) by integrating over the space coordinates 

p2tt poo 

N*(E) = dip rdrN*(r,E). (3.9) 
Jo Jo 

Since we have the integral representation (j37Sj) for Ag v (r, Q), it is straightforward to calculate directly the perturbation 
of DOS, AN^(E) = N V (E) - V 2D N^ induced by the Aharonov-Bohm potential. Here iVf = M/{2irh 2 ) is a free DOS 
of 2DEG per unit area and Vid is the 2D volume (area) of the system. Firstly integrating Ag v (r, Q) over space, one 
obtains^ 

2T , f°° , a , «x sin™ [°° dy f°° , cosh(2ry - l)x 

d<p rdrAg v (r,Q) = —+ — iU / dx V / ; . (3.10) 



o Jo 



Q 2 Jo cosh y Jo cosh x 



Then integrating over x and y we obtain that 



27T POO 



,1s I rdrAg v (r,Q) = 7rvil Q2 v) . (3.11) 



o Jo 
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Returning to the real q axis we reproduce the usual Aharonov-Bohm depletion of the DOS with respect to the free 
DOS V 2D N^: 

ANS(E) = -~\r,\(l-\ V \)5(E). (3.12) 
Writing Eq. fl3.12[) we have also included a case of the opposite field direction. 



C. The local density of states 



Now we come back to the LDOS p.6[) . As in the case of the full DOS, it is convenient to consider the excess LDOS, 
AN$(r, E) = N*>(r, E) - ivf. Then the value AN${r, E) can be obtained by calculating the function Ag v (r, Q) given 
by Eq. (|3.8[) and substituting the result to Eq. (13 4[) . The analytic continuation Q — > — iq described by Eq. (|3.6[) has 
to be done at the very last step of the calculation. 

Thus our purpose is to derive a simple representation for the function Ag v . First, one can rewrite it in the form 



Ag v (v ) Q) = - I 
Jo. 



dw 



dAg v (r,w) 
dw 



where we used that g v (r,oo) = 0. Differentiating Q3.8P we get the integrand of the last expression 



dxcosh[(277 - l)x] / dy coshye' 2Qrcoshxcoshy . 



dAg v (r,Q) 4 sin 7777 
dQ ~ ?r 1 

Using the integral representation for the modified Bessel function K v (x)^ 

/>oo 

K v (x)= dte- xcosht cosh ut 
Jo 

and the formula (2.16.13.2) from2£ 

J dxcoshbxK v (ccoshx) = ^K {l , +b)/2 K {u _ b)/2 Q 



we come to the equation 



dAg v (r, Q) 2sin7r77 



dQ 7T 

Integrating the last expression with Mathematica, we get 



rK n {Qr)K^ v (Qr). 



Ag v = 



sin irr/ 
8tt 



^(Qr) 2 - 2r >T 2 ( v - 1) 2 F 3 (1 - 77, 3/2 - 77; 2 - 277, 2 - 77, 2 - 77; Q 2 r 2 ) 



4 1 -"(Qr) 2 T 2 (-'7) 2*3 fa, 1/2 + r?; 2r?, 1 + 77, 1 + 77; Q 2 r 2 ) 



In 



(3.13) 



(3.14) 



(3.15) 



(3.16) 



(3.17) 



(3.18) 



- 4 ^ ) 3^(1, 1, 3/2; 2, 2, 2 - 77, 1 + 77; Q 2 r 2 ) + 1^(1 - v) + 

where p F q (ai, . . . a p ; b\, . . . b q ; z) is the generalized hypergeometric function and il>{z) is the logarithmic derivative of 
the gamma function T(z). After the analytic continuation Q — > —iq is made only two terms in the square brackets 
(with the hypergeometric function itself remaining real) and logarithmic term contribute in ImAg^r, — iq), so that 



AN^(r,E) = A^ s {sin 2 (7r77)[f (77,97-) +F(l-T),qr)] - 1} 



with 



F(f],qr) = 



4"- 1 (< ? r) 2 - 2r T 2 (?7- 1) 



2 F 3 ( 1 - 77, - - 77; 2 - 277, 2 - 77, 2 - 77; -{qrf 



(3.19) 



(3.20) 
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Note that the last term — Nq in Eq. (|3.19l) arises from the logarithmic term of (|3.18[) . In the important limiting cases 
the last expression is greatly simplified. For example, in the limit qr 3> 1 one obtains 



AN*(r,E) = -Ni 



, sin(7rry) cos(2qr) 



qr 



(3.21) 



Also in the physically important case r\ = 1/2 the LDOS is expressed in terms of the sine integral, Si(x) = J Q X dtsint/t 
as follows 



AN^ /2 (r,E) = N^ 



;Si(2 9 r) - 1 



(3.22) 



Using the asymptotic of the integral sine22, Si(x) « ir/2 — cosx/x — sinx/x 2 + 0(l/x 3 ) for x ^> 1 from p. 221) we 
recover the previous expression (|3.2ip valid for r\ = 1/2 and qr 3> 1. In the opposite limit Si(x) = x for x <C 1, we see 
that ANf /2 (r,E) = N$[Aqr/-n - 1]. 

In Fig. [T] we show the dependence (|3.22p of the induced LDOS AA^ 2 (r,_E) on the distance from the center of the 
vortex r. We observe that in the case of nonrelativistic 2DEG the presence of the vortex induced the depletion of the 
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FIG. 1: (Color online) The normalized LDOS function AN^{qr)/N^ as a function of the dimensionless variable qr for r\ — 1/2. 

LDOS for small qr <C 1. The function ANf^ 2 (qr) crosses zero near qr w 1 and for qr w 1.5 the function it reaches the 

maximal value ~ 0.2Aq . In Fig. [2] we model a situation when the STM tip is positioned at some distance from the 
center of vortex and a bias voltage is swept to explore the energy dependence of the LDOS. To take into account the 
presence of a finite carrier density in 2DEG, we introduce a finite Fermi energy so that the LDOS at zero energy, 
£ = 0, corresponds its value at the Fermi level, i.e. q = \/2ME/h — > q = ^/2I\I(£ + /j)/H. To choose the appropriate 
units we set the distance scale r to be the order of the lattice constant. Then the energy scale, E = h 2 / (2Mr$) is 
the order of the bandwidth. The dimensionless variable qr can now be rewritten as follows qr = y/£/Eo + n/Eor/ro- 
The dependence of AN^(£) is shown in Fig. [2] for three values of r/r^: for r/ro = 0.5 - solid (blue) curve, for r/ro = 1 
- long-dashed (red) curve and for r/ro = 10 - short-dashed (black) curve. The chemical potential is taken fi = Eq. 
We observe that only for the smallest ratio r/ro = 0.5 the values of the LDOS are significantly depleted below the 
free LDOS 7V S . As we saw, the depletion of the LDOS occurs for qr < 0.5. Since the presence of the Fermi surface 
makes the value of q large, so that the region of small qr is accessible only for r <^ rQ. Indeed, we observe that only for 
the smallest ratio r/ro = 0.5 the values of the LDOS are depleted to a half of the value of free LDOS Nq . Since the 
realistic values of the vortex core size are at least of the order of magnitude larger than ro, this implies that the region 
of a significant depletion of the LDOS is not accessible experimentally. Still, due to the slow decay of AA^ ~ 1/r 
even for r/ro — 10 the amplitude of AN$ oscillations is of order of 0.05Aq, so that this behavior can be probably 
observed experimentally. 



IV. RELATIVISTIC CASE 



In this section we examine the density of states for Dirac particles in the potential of the infinitesimally thin 
solenoid. To avoid formal complications we consider the physical regularization of the problem with the magnetic 
field concentrated in a thin cylindrical shell of small but finite radius R and take the limit R — > at the end of 
the calculation. All answers are presented in the form convenient for comparison with the case of the Schrodinger 
equation. 
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FIG. 2: (Color online) The normalized LDOS function AN%(£)/N§ as a function of energy £ for three values of r/ro = 0.5, 1, 10 
and /J, — Eq. 



A. Solutions of the Dirac equation in Aharonov-Bohm potential and general representation for LDOS 

The Dirac equation (I2.7[) in the field of the regularized vortex (I2.10[) , (|2.1ip (see also (lA5p ) is solved in Appendix lAl 
In our consideration we follow Refs.|21.22. The profile (| A5|) implies that for r < R the particle obeys the Dirac equation 
(|A7|) for a free particle, while for r > R the particle moves in the field of Aharonov-Bohm vortex (|A8I) . Accordingly, 
for r < R the squared Dirac equation (|2.8p is equivalent to free Schrodinger equations for the components of the 
spinor 

_ ( Mr) 



Notice that in the Appendix [A] the definition (|A2|) for ^2 explicitly includes the factor i. For r > R the components 
'l'(r) satisfy the Schrodinger equations with Aharonov-Bohm potential and the commutator (|2.9p is singular at r = R: 

i[D 1 ,D 2 ] = -l5(r-R). (4.2) 

The solution of the problem can found by matching the solutions obtained in the domains r < R and r > R [see e.g. 
(jAlip ], The radial components of the spinor \I/(r) have to be continuous: 

^iCR + o) = iM#-o), MR + o) = MR-Q), (4-3) 

and the singularity of the commutator (|4.2p is taken into account by a condition on the derivatives 

i>' 1 (R + o)-^' 1 (R-o) = ^MR), i/>' 2 (R + o)-il>' 3 (R-o) = -^MR)' ( 4 - 4 ) 

We stress that in contrast to the Dirac equation case, for the nonrelativistic case when the solution (|3.ip is obtained 
using the same regularization procedure, both the wave function and its derivative should be continuous. This as we 
saw from Eqs. (|2.8p . (|2.9p and (|A3I) is related to the pseudo-Zeeman term. 

After the limit R — > is taken we obtain for the case £ = 1 the following solutions: 



m{) VinE(k)\ ±ie i ^^E(k)-AJ lm+v] (kr) J' [ ^ 



for positive value of the energy E = E{k) = y/ (hvFk) 2 + A 2 and m^0; 



_ / fc / e*"- *> V^(fc) - A J |m+i7 _ X | (fcr) \ 
mU V 47r£; ( fc )V Tie^V^) + AJ |m+ ,|(fcr) ' 1 ^ 
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for negative value of energy E = —E(k) and m 7^ 0. Here the upper and lower sign of the second spinor's component 
corresponds to m > and m < solutions, respectively. The m = solution turns out to be special, while the upper 
component is regular at r = 0, the lower component diverges as J- V (kr) ~ r~'': 





f-AJ] 




- AJ_ 




- A J; 



for £ = £(fc) and 

for E = —E{k). One can check this property by using the matching conditions for ipi(r) component and then finding 
■02 (j) from Eq. (|A8b[) . On the other hand, one can discover this singularity in the limit R — > by directly analyzing 
the matching conditions for ip2{f) and then finding nonsingular 4'i( r ) from Eq. (|A8a[l . Comparing Eqs. (|4.7[) and (|4.8[) 
we observe that a singular at r — zero-mode solution E(k) — A is hole-like, because a singular electron-like solution 
vanishes due to the \/ E(k) — A factor. Concluding the discussion of the solutions for the case £ = 1 we stress that 
for the opposite field direction it is the upper component of the spinor which is singular in the solution equivalent to 
Eqs. (|4.7j) and ()4.8[) . Moreover, the nonvanishing zero mode is now electron- like. As it was firstly noticed in Ref. Hi], 
this behavior under the change in the field direction breaks the symmetry under $ — > <f> + $0 (see also Ref. [3l] for a 
discussion). 

The set of solutions for the case C = — 1 is the following: 



# fr] -J k ( e-yg (fc) + AJ |m +t?| (fcr) \ 



for E = E(k) and 



y(-) (r) J k ( e-yg (fc)-AJ| m +t? |(fcr) \ 
m U V 47r£(fc) ^ ±^(™-DV^(fc)+AJ|„ l+I) _ 1 |(fcr) ; ' 1 • UJ 

for i£ = —E(k). The prescription for the upper and lower sign of the second spinor's component is the same as for 
C = 1- The m = solutions in the case £ = — 1 are the following: 



*<+>M = / k ( VE (k) + AJ_ v (kr) . 
l J V 47r£(fc) ^ ie-^y/E\k) ~ AJi_„(fcr) ' ' ' 



for E = E(k) and 



(-) / k f y/E( k) - AJ-^ (fcr) 

W V 47r£(fc) V -ie-^s/Wk) + AJ^kr) 



(4.12) 



for i? = —E(k). We observe that in this case it is the upper component of the spinor which is singular at r = and 
nonvanishing is the electron-like zero mode. 

Now we are at the position to construct the GF for the Dirac fermions using the presented above solutions. However, 
before going to this, we should stress that as shown i n 21 ' 22 the field configuration with B confined to the surface of 
a cylinder of radius R is not essential for the main result. In simple words, an more general and not singular at 
the origin potential can be considered as a set of concentric shells. Or putting more formally, a profile function h(r) 
should satisfy the condition given below Eq. (|2 . 1 1 1) which excludes delta-funtion at the origin. 

The eigenfunction expansion for the retarded Dirac GF's now includes both positive and negative energy solutions 

G v (r,r,E + l O ) Q-J o E _ E{k) + tQ + E + E{k) + M J' ^ 

where depending on the sign of £ the sum is taken over either £ = 1 or £ = — 1 solutions. Accordingly, for the diagonal 
matric elements of G® , for instance, in the £ = 1 case we obtain 

E + A f°° kdk 



G° 11 {vy,E + iQ:X=l)= E ( ~t% I 



q 2 — k 2 + iOsgnE 

(4.14) 

m=—oo 



00 



9 



and 



E - A f°° kdk 



OO 

J- e m ^-^J lm+nl (kr)J lm+1ll (kr') (4.15) 



X 

7TL— — OO 

E - A f 00 fcdfc 



; [J_^(fcr-)J^(A;r') - J^fcr) J^kr')], 



2i:(hv F ) 2 J q 2 - k 2 + iOsgn£ 
where q 2 = (E 2 — A 2 )/ (hvp) 2 . Then the GF with coinciding arguments, r = r' acquires the following form 

E -\- A 

G° u (r,r,£ + iO;C=l) = , s 2 fln(r,g), 
G° 22 (r,r,^ + iO;C=l) = 2 f (fa; ^ )2 K(r, g) + /,(r, g)] 



(4.16) 



for C = 1) an d 



G° u (r, r, £ + zO; C = -1) = ^ [ff,(r, <?) + /„(r, <?)] 

A 

G® 22 (r,r,E + iO;C= -1) = 2 ^ hvp y 9n( r > g)» 



(4.17) 



for £ = — 1. Here the function g v (r,q) is related to the function ^(r, Q) defined in Eq. (|3.5p by analytic continuation 
<7 + iOsgnE' — >• z = iQsgnE. Note that while the function <7„(r, Q) is identical for both nonrelativistic and relativistic 
cases, the real momentum q function g n (r,q) has a different analytical properties in these cases reflecting the fact 
that in contrast to the nonrelativistic case the relativistic spectrum contains positive and negative energy branches. 
Similarly, the function 

r°° kdk 

f v (r,Q) = -J Q2^i J -v( kr ) ~ J n( kr )\- ( 4 - 18 ) 

is the analytic continuation to the imaginary axis of the function /^(r, q). The f v contribution to G D originates from 
the zero mode solutions of the Dirac equation. Using the integral [see Eq. (2.12.32.12) from^] 

00 kdk 

J 2 v {kr) = I v {cr)K v {cr), (4.19) 



o c 2 + k 2 

where I v (cr) is the modified Bessel function and K v {cr) is the Macdonald function, we obtain the following simple 
result for f n 

Mr,Q) = - 2 -^K*(Qr). (4.20) 
B. The density of states 

Due to the matrix structure of the GF and the presence of the valley degree of freedom, £ = ±1, in contrast to the 
nonrelativistic expression (I3.9[) . the full DOS in the relativistic case involves not only the integration over area, but 
also summing over the diagonal components of the GF and K± valleys. For a better understanding of the final result 
for the full DOS, it is instructive to consider separate expressions for the DOS corresponding separate K± points 



2 p27T />00 

p%(E,C) = / dip rdrlm [trG?(r,r,£ + »0;0] ■ (4-21) 

n Jo Jo 



Recall that for 77 = the free DOS of the Dirac quasiparticles per spin and one valley is equal to po(E) — \E\9(E 2 / A 2 - 
l)/(2nfi 2 v 2 F ). 
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Using Eqs. (14.16[) and (|4.17[) one can express the DOS (I4.2ip via the functions g v (Y, 2) and /^(r, Q) defined by 
Eqs. p.5[) and ([4.181) . respectively, as follows 



,2* f J (4-22) 

x dip rdrIm[2Eg 7l (r,Q^-iqsgnE + 0) + (E-CA)fr l {r,Q^-iqsgnE + 0)]. 
Jo Jo 

The integral of used in the nonrelativistic case function g^r, Q), or to be precise of the function Ag, v over the 
space coordinates is calculated in Eq. (|3 . 1 1 [) . The corresponding space integration of the function f n (r, Q) given by 
Eq. (|4.20p produces the result 

dip j rdrf v (r,Q) = —gf. (4.23) 

Notice that since f v = for r\ = 0, there is no need to introduce a function Af, r Having these results for the spatial 
integration we can calculate the imaginary part. Since both Eqs. (|3 . 1 1 [) and (|4.23|) depend on Q 2 , the imaginary part 
is obtained by using a simple prescription Q 2 — > — q 2 = —[(E + i0) 2 — A 2 }/(Hvf) 2 ■ This gives the final expression for 
the perturbed DOS, Ap%(E, () = p®(E, Q - V 2D Po(E) by the Aharonov-Bohm vortex 

Ap»(E, = -\\v\(l - \r}\)[6(E - A) + S(E + A)] + \ V \5(E + Csgn(r?)A). (4.24) 

The first term of (I4.24[) which is ~ — |tj|(1 — |f?|) originates from g v part of Eq. (|4.22[) and resembles the nonrelativistic 
result (|3.12j) . For A = it turns out to be twice larger than (|3.12j) simply because in Eq. (|4.24|) we summed over 
the diagonal components of the GF GP which are related to the two sublattices of graphene. The last term of 
(|4.24p originates from f v part of Eq. (I4.22|) and thus is related to the zero mode solution of the Dirac equation. The 
mentioned above fact that the singular component of the zero mode solution is hole-like for £ = 1 [see Eq. (|4.8p ] and 
it is electron-like for £ = — 1 [see Eq. (|4.11[) ] finds its reflection in the asymmetric form of the last term of Eq. (|4.24p 
which also corresponds to the holes (electrons) for £ = 1 (£ = — 1). In accordance with the behavior of the zero 
mode solutions described in Sec. 11V A[ when the direction of the field is reversed the expressions p v (E,( = 1) and 
p v (E, ( = —1) arc interchanged. The full excess DOS 

AN»(E) = Ap^(E, C = 1) + Ap°(E, C = -1) = V 2 [S(E - A) + 5(E + A)} (4.25) 

is obviously symmetric in energy. In contrast to the nonrelativistic case, in the Dirac case the Aharonov-Bohm vortex 
induces the excess of the states which is related to the presence of the last term of Eq. (|4.24|l and caused by the zero 
modes. We note that in Ref. [l3| the corresponding term of Eq. (|4.24p has a wrong sign. The positiveness of AN^(E) 
can also be understood by the following simple argument. For A = the free DOS p$(E = 0) = 0. Therefore, since 
the DOS has to be positive, the value AN^(E) should also be positive. 

C. The local density of states 

Now we investigate the LDOS for the Dirac case. While it was useful to consider each valley separately, especially 
because field theoretical studies of the problem often involve only one unitary inequivalent representation of 2 x 2 
gamma matrices [see e.gJ^], the LDOS measurement picks up both valleys together. On the other hand, LDOS 
distinguishes sublattices. Thus we consider separately the LDOS for A and B sublattices which are defined as follows 

(r, E) = -~Im [G vll {r, r,E + iO;(=l) + G, ;11 (r, r, E + zO; C = -1)] , 

I (4.26) 
N^ B \r, E) = Im [G^r, r, E + i0; C = 1) + G, )22 (r, r, E + i0; ( = -1)] . 

Again we consider the perturbation of the LDOS induced by the Aharonov-Bohm potential, AN^ (r, E) — 
N^ (A ' B) (r,E) - A D(j4 ' S) (r,£). Using Eqs. (|4T6| and (j4T7| we rewrite the LDOS's AN^ (A - B \r, E) in terms of the 
functions Ag v (r, Q) and /,,(r, Q) as follows 

AN^ A > B \r, E) = ~ 2 f ^ j 2 Im[2Aff n (r, S -> -iqsgpE + 0) + / n (r, Q -> -iqsgnE + 0)]. (4.27) 
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Here the upper (lower) sign corresponds to A (B) sublattice. To obtain the final expression for LDOS we should make 
the analytic continuation Q — > —iqsgnE. For the function Ag v (r, Q) given by Eq. (|3.18p one obtains 

ImA 5t) (r, Q -> -iqagnE) = ~ ^|— {sin 2 7nj[F(jj, qr) + F(l - 77, qr)} - 1} , (4.28) 

where the function F(r],qr) is defined in Eq. (|3.20p . One can easily see that because in the nonrelativistic case the 
analytic continuation to the real momentum Q — > —iq is different from the analytic continuation Q — > —iqsgnE in 
the relativistic case, the function (|4.28l) differs from the function (|3.19[) . Using the relationships^ 

K v {ze^) = -i^e-^HW(z), H^(z) = J v {z) - iY v {z) (4.29) 

between the Macdonald function of the imaginary argument, the Hankel function of the second kind Hi (z) and the 
Bessel function of the first J v (z) and second Y v (z) kinds, we obtain that 

K 2 (±iz) = ^ e ^[Y 2 (z) - J 2 (z) ± 2iJ n (z)Y r ,(z)}, (4.30) 

where we used the property K v (z*) — K*(z). Accordingly, the analytic continuation of the function /^(r, Q) from 
Eq. (|4~20| takes the form 

Im/„(r, Q -> -iqsgnE) = - - sm7T ^ s S nE { sinjTri [ Y 2 (qr) - J*(qr)} - 2cos7rr ? J,(«r)Y I ,( 9 r)}. (4-31) 
Substituting (|4.28|) and f|4.31|) in Eq. (|4.27j) we arrive at the final main result 



AiV^ fl )(r,^) = iV D ( J E) (l±§)» 




sin 2 (7r77)[F(?7, qr) + F(l — r], qr)] — 1 



-^K,V) - Jfar)] ~ ^Mqr)Y v{ qr)\ q = 2^=^ 



(4.32) 



where Nq > (E) — \E\/ (2irti 2 Vp) is free DOS of the Dirac quasiparticles per spin and one sublattice (or valley) for 
A = 0. The first part of Eq. (|4.32l) which includes F and —1 is identical to the nonrelativistic expression (|3.19[) . while 
the second part of Eq. (|4.32p with Bessel functions originates from the zero mode contribution. 
In the limit qr 3> 1 the last expression acquires a simple form 

A<<-», r , E , . (l ± §) , (f! - l) '^y . (4.33) 

Comparing Eqs. (|4.33p and p.2ip we conclude that in the Dirac case the impact of the vortex is more localized than 
in the nonrelativistic case. 

As we saw in Sec. IIII Cl in the physically important case rj = 1/2 the expression (|3.19p is significantly simplified to 
the result (|3.22p . The same remains true for Eq. (|4.32l) . because half-integer Bessel functions are expressed in terms 
of the elementary functions and we obtain that 



AN^- B \r,E)=N»(E)(l±^) 0(^-1 



*Si(2qr) 1 + ^» 
7r ' irqr 



(4.34) 



From Eq. (|4.34l) we immediately observe the main difference between the relativistic and nonrelativistic cases. The 

presence of zero modes causes a positive divergence of the LDOS, AN^ A ' B \r, E) ~ l/qr for qr <C 1 near the center 
of the vortex. Integrally this results in the excess of the states in the full DOS (|4.25p . Using the asymptotic expansion 
of Si(x) given below Eq. (|3.22p we recover the previous expression (|4.33p valid for rj = 1/2 and qr ^> 1. 

In Fig. [3] we show the dependence (|4.34l) of the induced LDOS ANy% ' B \t, E) on the distance from the center of 
the vortex r for A = 0. We observe the features expected from the analytic expressions such as the excess of the 
LDOS for small qr < 1 and faster than in 2DEG decay of ANR 2 — 1/r 2 for qr > 1. 

In Fig. 2] we consider a situation similar to Fig. [21 We fix the distance at r — 10r and plot the energy dependence 
of the relativistic LDOS (|4.34p . We consider the most interesting case of undoped graphene with zero carrier density. 
In contrast to the 2DEG, in graphene is easily tuned to this regime. Again, we introduce the distance scale rg of 




FIG. 4: (Color online) The normalized LDOS function AN° (A ' B) (E) /N$ (E) as a function of energy E for r/r = 10, A = 0.1E 
and r\ = 1/2. 



the order of the lattice constant. Then for the Dirac case the energy scale is E = hvp/r , and accordingly the 
dimensionless variable is qr — ^/(E/Eq) 2 — (A/ Eo) 2 r /rrj. To make sublattices inequivalent we introduce a finite gap 
A = O.lE'o which introduces the asymmetry between the the LDOS on A and B sublattices. The LDOS AN^^{E) 
shown as a solid (blue) curve has a sharp peak near E = A and practically no peak at E = A, while the LDOS 
AN^I~^\e) shown as dashed (red) curve has a peak near E = — A and no peak at E = A. When E is increasing the 
LDOS very quickly reduces to its value for the Dirac system in the absence of the vortex. We note that this behavior 
of the LDOS is seen for the large value of r/r = 10, i.e. far from the center of the vortex, indicates that it should be 
possible to observe this excess of the LDOS in experiments on graphene. 

Finishing our discussion of the Dirac case, we mention a recent work22, where the effect of the vacuum polarization 
in the field of an infinitesimally thin solenoid at the distances much larger than the radius of solenoid was studied. 
Constructing the GF the authors neglected the delta-function (14.41) motivated by the fact that they are interested in 
the regime r 3> R. One of their interesting conclusions is that for r\ = 1/2 the induced current is zero. 



V. CONCLUSIONS 



The main motivation of this work was to address a natural question whether one can distinguish graphene from 
2DEG measuring the LDOS near the Abrikosov vortex penetrating them. Just by comparing Figs. [2] and [4] we find a 
positive answer. Indeed, by putting the STM tip close to the vortex, for example, at the distance, r ~ 10ro, where ro 
is the lattice constant, we observe that in 2DEG the value of the LDOS will be close to its background value Nq in the 
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free system, while in undoped graphene one should observe a strong LDOS enhancement near the Fermi level. Or to 
be more precise, if there is no gap A in the quasiparticle spectrum, the peak of the LDOS should be observed at the 
Fermi level, or if A 7^ the peaks should be observed at the energies E — ±A and the sign of the peak energy should 
depend on the sublattice. This peaked behavior of the LDOS with AN^ A ' B \r, E) ~ l/qr for qr -c 1 reflects the 
specific feature of the Dirac fermions such as the presence of the divergent as 1/ \fr at the origin zero mode solution 
of the Dirac equation. Thus the observation of this feature in STS measurements would contribute to the expanding 
list of the experimental manifestations of the Dirac fermions in graphene. Our second conclusion is that while in the 
nonrelativistic case the presence of the Aharonov-Bohm vortex leads to a depletion of the full DOS, in the Dirac case 
the full DOS is enhanced. This result can likely be checked by analyzing the STS maps. We remind at the end that 
in this paper we did not consider the effect of disorder in the presence of magnetic field which will definitely affect 
the behavior of the LDOS in magnetic field^ and should be taken into account in the analysis of the STS data. 
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Appendix A: Solution of the Dirac equation 

A positive energy solution of the time-dependent Dirac equation has a form r) = cxp(—iEt/h)^f(r) with a 
two-component spinor \&(r) satisfying the time- independent Dirac equation (|2.7[) : 



ihv F o x \ dx + i Y c A ^) + iKv F C,a 2 [d 2 + i-^A" 1 ) ~ Ac73 \ *( r ) = °- ( A1 ) 

where £ = ±1 distinguishes two unitary inequivalent representations of 2 x 2 gamma matrices. It is convenient to 
denote the components of two-component spinor 'J'(r) 

•<"-(&$) < A2 > 

with the factor i explicitly included in the definition of the lower component. Then we can rewrite Dirac equation 
(|Aip in the components as follows: 

(E - A)Vi(r) - HvpiDx - iCAOlkW =0, 

hv F {D 1 +iCD 2 )i; 1 {v) + {E + A)ij 2 {v)=Q. ( ' 

Since we consider a cylindrically symmetric configuration of the field with a vector potential A = e ip A tp {r) 1 the system 
13)) has to be rewritten in the polar coordinates (r, ip): 



(A4) 



As discussed in Sec. [TT1 to analyze the problem with a singular at r = Aharonov-Bohm potential (|2 .4[) one has to 
do a self-adjoint extension of the Dirac operator, see e.g. Refs. [T2IIT4II27I . To avoid this complication we consider a 
regularized field configuration (|2.10|) suggested in Refs. l2lU22i with the profile function (|2.1ip . so that 

eAyjr) = f0, r < R, 
he I rj/r, r > R. 
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From now on, we consider the specific case ( = 1 and seek for a solution of Eq. (IA4[) in the following form 



Mr) = e* (m -^Mr), M') = e mip Mr) 
for r < R we obtain a system of the radial equations for a free Dirac particle 



(E - A)ipi(r) - hv F 
to — 1 



d to 
dr r 



d_ 

dr 



Mr) =0, 
Mr) + (E + A)Mr) =0, 



while for r > R we have 



(E - A)M r ) - ftvp 
d 1 



dr 



-(m + r) — 1) 



d 1 1 

— + -{m + r]) M r ) =°> 

ar r 

^i(r) + (£ + A)V> 2 (r) =0. 



(A6) 

(A7a) 
(A7b) 

(A8a) 
(A8b) 



One can obtain from the systems (|A7[) and (IA8|) that the spinor components satisfy the following second order 
differential equations: 



or 



1 d 



r dr 
1 d 
r c?r 



M r ) 



lr' 
d 2 

-^2^2{r) + z—^ r ) 



E 2 



(m-1) 2 

r 2 (Hvf 
m 2 E 2 - A 2 



Mr)=0, 



Mr)=0, 



for r < R and 



d 2 , . Id,,, 
ar r ar 



(owf) 2 
(to + ?7 - l) 2 E 12 - A 2 



-r^Mr) + --rMn - 

dr z r dr 



(to + riY 



(hv F ) 2 

E 2 - A 2 
(hv F ) 2 



Mr)=0, 
Mr) = 



(A9a) 
(A9b) 

(AlOa) 
(AlOb) 



for r > R. The solutions of Eqs. (|A9[) and (|A10[) are expressed in terms of the usual Bessel functions. For example, 
for the solutions of the equations (|A9a[) and (jAlOap for the component ifji ( r ) are given by 



M r ) = C mJ\m-l\(kr), r < R, 

M r ) = A mJ\m+r 1 -l\{kr) + B m J_| m+ ,_1| (Ar), 



r > R, 



(Alia) 
(Allb) 



where A m , B m , and C m are constants and the J's are the Bessel functions. The solution (|Alla[) is standard due to 
normalizability and the absence of delta function at r — 0. The coefficients A m , B m , and C m to be found from the 
matching conditions (|4.3[) and (14.41) . To find the second component, one can substitute the result for ipi(r) in 

Eq. (|A8b[) . Or, equivalently, one can start from Eqs. (|A9bl) and (|A10b[) for the component ip2(r) which also have a 
solution in the form (|A11|) . find the corresponding constants from the matching conditions and then use Eq. (|A8al) to 
obtain ipi (r) . Finally, the overall factor before the solution is determined by the normalization condition 



2tt 



dip / rdr^ m , (r, ip; k')V m (r, <p; k) = S(k - k')S„ 



(A12) 
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Notice that the angular momentum operator does not commute with the Dirac Hamiltonian. 

In what follows we also consider the behavior of the results under the reversal of the field direction which corresponds to 
the negative values of rj. 



